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Abstract 


The resolution requirements for direct numerical simulations of transition to turbulence are 
investigated. A reliable resolution criterion is determined from the results of several detailed 
simulations of channel and boundary-layer transition. 
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Introduction 


In a numerical simulation of fluid dynamics (as in an experiment) it is surely far less 
difficult to report what is observed than it is to judge the extent to which the observations 
are meaningful. Most scientists have been inculcated with a healthy skepticism about the 
precision of experimental data, but perhaps because of their comparative novelty (and, yes, 
perhaps also because of the increasingly seductive allure of computer graphics) numerical 
results receive less critical acceptance. 

The conventional wisdom on the failings of coarse-grid numerical simulations is that 
artificial viscosity produces excessive smearing of the flow field. However, in certain cir- 
cumstances the real flaw of coarse-grid computations is just the opposite - inadequate 
numerical resolution induces spurious oscillations in the flow. In direct numerical sim- 
ulations of transition from laminar to turbulent flow, this latter effect is the principal 
one. 

Transition problems can be divided into two broad categories: bifurcation problems, 
in which a regular sequence of instabilities arises as the Reynolds number is varied, until 
at large enough Reynolds number the flow is chaotic; and abrupt transition problems, in 
which a flow which is unstable at a fixed Reynolds number passes through a rapid sequence 
of instabilities leading inevitably to turbulence. Three-dimensional, nonlinear numerical 
calculations of the time-dependent Navier-Stokes equations have been used extensively in 
the last 5 to 10 years to investigate many important features of transition from laminar to 
turbulent flow. Taylor-Couette (Marcus (1984); Streett and Hussaini (1987)) and Rayleigh- 
Benard (Curry et al (1984)) flow have been the focal points of numerical simulations of 
bifurcation problems, whereas free shear layers (Metcalfe et al (1987); Lowery et al (1987)) 
and channel and boundary-layer flows (see the review by Zang and Hussaini (1987)) have 
been the prime targets of abrupt transition simulations. 

Orszag and Kells (1980) have characterized the artificial oscillations that arise in poorly 
resolved flows as “numerical turbulence”. In bifurcation problems numerical turbulence 
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may lead one to conclude erroneously that a given flow is chaotic (Curry et al (1984)). In 
abrupt transition problems, it could lead one to mislocate the instant and region in which 
turbulence originates as well as to mis-identify the precise instability which induces the 
turbulence. 

Consequently, numerical simulations must be conducted with extreme care, with special 
attention paid to providing sufficient numerical resolution, i.e., enough grid-points. The 
resolution requirements for abrupt transition problems are particularly stringent, since as 
an unstable flow evolves towards a turbulent state the requisite resolution demands increase 
dramatically. In this paper we illustrate this point in some detail for some representative 
simulations of transition to turbulence in channel and boundary-layer flow. Moreover, 
we furnish some guidance on the level of accuracy that is necessary to obtain meaningful 
results. 


Numerical Approximation 


The channel flow geometry is illustrated in Figure 1. Lengths are scaled by the channel 
half-width h , velocities by the centerline velocity u c , and the density by its constant value. 
The Reynolds number Re is based on these scales, i.e., Re = u c h/u, where u is the 
kinematic viscosity. The flow is assumed to be fully developed with the parabolic mean 
flow profile uo(y) = 1 — y 2 . 

In dimensionless variables, the incompressible Navier-Stokes equations in rotation form 
may be written 

^ +w x u + V(^|u| 2 ) + Vp = i/Au. (1) 

V • u = 0 (2) 

where u = (u, v, w) is the velocity, w = V x u is the vorticity, and p is the pressure. The 
independent variables are assumed to be periodic in x and z , with no-slip conditions on 
u at the walls y = ±1. 

Two related numerical methods have been used to obtain the results reported herein. 
Both use a Fourier collocation technique for the spatial discretization in the x and z direc- 
tions. One uses a second-order finite-difference scheme in y, and the other uses Chebyshev 
collocation for that direction. A Fourier-finite-difference method (FFDM) was devised by 
Wray and Hussaini (1984) for direct simulation of boundary-layer transition and by Moin 
and Kim (1982) for large-eddy simulation of channel flow; the latter version was used by 
Biringen (1984, 1987) for his simulations of channel flow transition and for some of the 
present results. The particular Fourier-Chebyshev method (FCBM) used here is the one 
developed by Zang and Hussaini (1986). The FCBM code uses the standard Chebyshev 
grid in y, 

y = -cos £, (3) 
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with the grid-points uniformly distributed in the computational variable £ e [0, 7 r|. The 
FFDM code uses the hyperbolic tangent transformation, 


y = —tanh{£ tank x a) 
a 


( 4 ) 


with the stretching parameter a = 0.94175 and £e[— 1, 1]. Both methods employ an explicit 
time discretization for the advection and horizontal diffusion terms and an implicit time 
discretization for the normal diffusion and the pressure gradient terms as well as for the 
continuity equation. The FCBM uses a splitting, or fractional-step, technique for the time 
advancement: in the first half-step the advection and diffusion terms are accounted for, 
and in the second half-step a pressure correction is applied so that the final velocity field 
is divergence free. Boundary conditions are applied at the intermediate step to all three 
velocity components and in the pressure correction to the normal velocity component. The 
FFDM is unsplit, i.e., all terms are accounted for simultaneously in a single step. Both 
codes employ a staggered grid for the pressure so that no artificial boundary conditions 
are required. A thorough discussion of the available spectral algorithms for channel flow 
simulations is provided in Chapter 7 of Canuto, et al (1988). 

Recently, Zang (1989) demonstrated that in numerical simulations utilizing collocation, 
it is much better to use the Navier-Stokes equations in the skew-symmetric form, given by 


du 

~dt 


+ • Vu + • (uu) + Vp = v Au 

I It 


( 5 ) 


than in the commonly used rotation form given by Eq. (1). Several numerical examples 
were presented, including one for channel flow, which demonstrated the superiority of the 
skew-symmetric form. This result is quite recent and post-dates some of the calculations 
reported in the present paper. Unless otherwise noted, the computations reported here 
utilized the skew-symmetric form. 

In these algorithms, both physical space and spectral space representations of the de- 
pendent variables are employed. The velocity has the Fourier series representation 


u(x,t) =u 0 (y)+ X) H n kxik ,{y,t)e^ ax ^ pz 1, (6) 

k x =-N x /2 k x =-N x / 2 

where u(y) = (uo(y),0,0) , and a and /? are the fundamental wavenumbers in the stream- 
wise and spanwise directions respectively. The fundamental wavelengths in these directions 
are given by L x = 2n/a and L z = 2tt//J, respectively. The imposed periodicity lengths 
are L x and L z . The integers k x and k z label the Fourier wavenumbers in the numerical 
representation with respect to the fundamental wavenumbers a and /?. The velocity also 
has the Chebyshev series representation 



n=0 


( 7 ) 
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For later reference we list here a useful measure of the strength of a given Fourier 
harmonic: 

Ek„k.{t) = /_ x I *) \*dy/E 0 , (8) 

where 

dk = 2 — 6k o. (9) 

This quantity is the kinetic energy of the mode normalized by the kinetic energy of the 
mean flow (E 0 ). The constant d* is used to account for the symmetries in the Fourier 
harmonics. Similarly, a measure of a given Chebyshev component is 

E n (t) = f f |u n (x,z,t)| J dxdz/(L*L,£ 0 ). (10) 

JO Jo 

Note that Eo,o(t) and Eo(t) measure the kinetic energy of the mean flow correction (see 
Eq. (6)). 

The initial conditions for the transition simulations reported in this work have the form 

u(x,0) = 3fce{uo(y) + c 2 r» e <# u 2 u(y)e ia * + ^e 3D [u£,(y) «<(«■+**) + u^y) (11) 

where u 2 £>(y) and (y) are the least stable linear modes for the given real wavenumbers a 
and /?. These eigenfunctions are commonly normalized so that their maximum streamwise 
amplitudes are 1, and so that, when used in the form given in Eq. (11), these maxima occur 
at x = 0 and z — 0 for 0 — 0. Thus, e 2 p and e 3 d measure the maxima of the streamwise 
fluctuations of the 2-D and oblique disturbances relative to the mean centerline velocity, 
and 0 measures the phase shift between the 2-D and the oblique waves. 


A Re = 1500 Channel Transition 


The particular channel flow transition studied by Biringen (1984, 1987) is used here as 
the primary example of the resolution demands of transition simulations. The Reynolds 
number is 1500. The horizontal wave numbers a and 0 are both unity. The tempo- 
ral frequencies of the least rapidly damped modes are u> 2D = 0.326299 — 0.028206* and 
u 3D = 0.401293 — 0.028230*, as determined by a Chebyshev tau approximation to the Orr- 
Sommerfeld equation (Orszag (1971)). The initial amplitudes of the waves were c 2 z? = 0.11 
and e 3 £> = 0.10, and the phase shift was 0 = -37°. The numerical computations are made 
in a frame moving with a streamwise velocity of 0.36. This is the average of the phase 
speeds of the two- and three-dimensional waves. The use of a moving frame improves the 
temporal accuracy of the calculation. All calculations presented in this section employed 
the rotation form. 

A straightforward test of the resolution in the normal direction is the ability of a code 
to reproduce the time evolution of a linear eigenfunction. The computed growth rates and 
their errors at t = 30 are given in Table 1 for several calculations. The initial condition 
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for each case consisted of the mean flow plus an oblique TS wave at 0.01% amplitude. 
Since Fourier collocation is used in the horizontal directions and the linear solution is 
a pure trigonometric function in x and z, a mere 4 points in each of these directions 
suffices for essentially perfect resolution. The time-step was chosen sufficiently small so 
that the errors were dominated by the normal discretization. The Chebyshev calculation 
is accurate to better than five significant digits. The finite-difference calculations exhibit 
the expected second-order accuracy. Moreover, the finite-difference calculations err on the 
side of excessive damping of the wave. (This same effect was found by Malik, Zang, and 
Hussaini (1985) at a higher Reynolds number.) 


Algorithm 

N v 

Calculated 

Error in w< 

Exact 


-0.028230 


FCBM 

64 

-0.028226 

-0.000004 

FFDM 

50 

-0.032030 

-0.003800 

FFDM 

100 

-0.029230 

-0.001000 

FFDM 

200 

-0.028482 

-0.000252 


Table 1: Effect of Normal Resolution on Growth Rate 


The error in the N v = 50 FFDM growth rate might appear to be small; after all, it 
is less than 1% when normalized by |wsp|. Nevertheless, it portends uncomfortably large 
errors in the full transition simulation. These manifest themselves quite clearly in the time- 
dependent vortical flow field. Since the vorticity is obtained from the primitive variables 
by differentiation, it is a sensitive indicator of errors in the velocities themselves. 

At early stages in the transition process there are two pronounced structures in the 
flow: a predominantly streamwise vortex known as a lambda vortex and a strong detached 
shear layer which lies above it. These features have been documented experimentally 
by Williams, Fasel, and Hama (1984) and numerically by Krist and Zang (1987), among 
others. The bottom row of Figure 2 shows the results in the lower half of the channel 
from a well-resolved, 64 s FCBM calculation at t = 15. The two quantities displayed 
are the vertical shear, du/dy , and the streamwise vorticity, u x . The former is shown in 
the spanwise plane of symmetry, z/L M = 1/2, which is the so-called “peak plane”. The 
detached shear layer is located in the upper left half of the frame. The latter is given for 
a streamwise plane near the maximum vertical shear, in this case at x/L* = 1/2. The 
lambda vortex is reflected by the uppermost counter-rotating features near the center of 
the frame. An thorough description of the flow field for this case can be found in Krist 
and Zang (1987). 

The results from a 32 x 50 x 32 FFDM calculation and from a 32 x 64 X 32 FCBM cal- 
culation are shown in the top two rows of Figure 2. (The FFDM grid was chosen to match 
that used by Biringen (1984).) The two fully spectral results are in close agreement. The 
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coaxser one does, however, exhibit some slight oscillations in the streamwise vorticity con- 
tours. The errors in this calculation stem primarily from inadequate spanwise resolution. 
Errors in the FFDM calculation are much more pronounced. The detached shear layer 
clearly lags behind its true position. Moreover, there are spurious streamwise vorticity 
regions located above the lambda vortex. The principal errors here are those due to the 
discretization in the normal direction; errors from the inadequate spanwise resolution are 
present at about the same level as in the coarse FCBM calculation. Krist and Zang (1987) 
compared the FFDM results shown here with FFDM results on a 32 x 200 x 32 grid. The 
phase lag of the detached shear layer is not nearly so severe, and the spurious streamwise 
vorticity regions are smaller. Nevertheless the results are still markedly inferior to those 
from the 32 x 64 x 32 FCBM. Bear in mind that Table 1 indicates that with 200 normal 
grid-points, the error of the FFDM on the primary waves is less than 0.1%. This fact 
underscores the extreme resolution demands posed by numerical simulations of transition 
to turbulence. 

After t = 15, the flow field evolves rapidly and the resolution demands increase sub- 
stantially. Figure 3 presents the results for the coarse-grid FFDM, the 64 3 FCBM, and a 
well-resolved FCBM at t = 18.75. (The well-resolved results here and in Figures 6 and 
7 were the product of a single simulation in which the grid was refined (using spectral 
interpolation) whenever further resolution was required.) Note that the detached shear 
layer has intensified as has the streamwise vorticity. Notice also that the streamwise vortex 
has pinched closer towards the peak plane and that it has narrowed its dimensions in the 
spanwise direction. Hence, the flow field contains sharper gradients and more grid-points 
are needed. In fact, the spanwise resolution demands are more than twice as great as they 
were at t = 15. The 64 s FCBM result is plainly inadequate. The streamwise vorticity 
contains an obvious grid-scale oscillation, and problems are evident in the vertical shear 
near x/L* = 3/4. The 32 x 50 x 32 FFDM result is afflicted with severe oscillations in 
the streamwise vorticity, and the quality of the vertical shear has also deteriorated. The 
effects of the phase lag are still present. 

The Fourier-Chebyshev coefficients of the computed solution are a good indicator of 
the accuracy of a fully spectral result. The resolution spectra of several cases are shown in 
Figure 4. The resolution in the normal direction is measured by the quantity E n (t) defined 
by Eq. (8) , whereas the resolutions in x and z are measured by 


EkA*) = max E kxik , (<) 

*jr 

(12) 

E k At) = maxE fciit ,(f), 

(13) 


respectively. Consider first the 64 s case at t = 15. There is a decay of at least 8 orders of 
magnitude from the low wavenumber end to the high wavenumber end of the resolution 
spectra in each coordinate direction. Since the spectra are quadratic in the perturbation 
velocities, this corresponds to a truncation error of better than 0.01% in the primitive 
variables themselves. The resolution spectra in x and z have a mild curl at their tails. In 
the case of the x-direction this is due to the round-off error of these 32-bit calculations. 
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The curl of the tail in the 2 -direction is at a higher level than can be accounted for by 
round-off error. It signals an imminent loss of resolution on this grid, a claim supported 
by the 64 s results at t = 18.75 shown in Figure 3. The resolution spectra at t = 18.75 for 
the well-resolved case also exhibit an 8 order of magnitude decay. In neither the t = 15, 
64 s nor the t = 18.75, 96 x 128 x 162 contour plots is there visible evidence of oscillations. 
Thus, in both of these well-resolved cases the lack of visible evidence of oscillations in the 
vorticity field corresponds to a truncation error of less than 0.01%. 

The bottom right example in Figure 4 is for the under-resolved FCBM case at t = 
18.75. The resolution spectra suggest that the x-direction is resolved the best, and the 
2 -direction the worst. This is consistent with the middle row of Figure 3, which indicates 
severe oscillations in the spanwise direction and few problems in the streamwise direction. 

One must be careful in how the data is plotted when looking for small-scale oscillations. 
This point is driven home by Figure 5. Here the results of a 64 s FFDM simulation at t 
= 18.75 are plotted in two ways. In the top row, all points in the flow field are used 
in the plotting. In the bottom row, every other point (in the horizontal directions) has 
been discarded before plotting. Notice that the flow field appears to be much smoother 
in the bottom case. The explanation is that the most severe oscillations are of the two- 
point, sawtooth variety and these are exactly eliminated by the coarsening procedure. 
Since simulations of transition require large data bases, this type of coarsening for data 
reduction may be adopted in the interests of conserving storage. Nevertheless, it can 
present a very misleading picture of the reliability of the results. In addition, the use of 
plotting routines which deliberately smooth the data should be avoided. 

Incidentally, a comparison of the top row of Figure 5 with the top two rows of Figure 
3 confirms that the spurious streamwise vorticity regions are indeed due to the FFDM’s 
inadequate normal resolution. 

The effects of inadequate resolution are even more pronounced at later times. Figure 
6 shows the computed flow at t = 22.5. The FFDM results exhibit only the barest re- 
semblance to even the large-scale vertical shear, and the streamwise vorticity is extremely 
corrupted. The coarse-grid FCBM flow contains large-amplitude grid-scale oscillations in 
the streamwise vorticity, but the vertical shear is still reasonable. The resolution spectra 
for the well-resolved case is given in the lower left frame of Figure 4. Once again the 8 
order of magnitude decay correlates well with the absence of visible flow-field oscillations. 

The t = 27.5 results are given in Figure 7. The instantaneous vorticity field from the 
FFDM is now completely unreliable. The large-scale features of the coarse-grid FCBM 
vertical shear are still correct. However, there is substantial contamination of the small- 
scale features in the bottom quarter of the plot, near the lower wall of the channel. The 
well-resolved case still has a smooth flow field and its resolution spectra (not shown) are 
again indicative of a 0.01% truncation error. 

From t = 25 onwards the flow is undergoing the final stages of transition to turbulence. 
As documented by Zang and Krist (1989), the principal flow field activity shifts markedly 
downward towards the near-wall region. It is precisely here that the coarse-grid results are 
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most unrealistic. Hence, they are effectively useless in addressing the physics of the final 
stages of transition. 

The preceding comparisons indicate that loss of resolution on the small-scale flow fea- 
tures leads eventually to loss of resolution on the large-scale features. Quantities such as 
the resolution spectra appear to be the best measure of the resolution for spectral meth- 
ods. Unfortunately, the resolution spectra are not quite as meaningful in a direction which 
is discretized via finite differences, since the spectra are not so directly related to the 
numerical approximation. 

The 0.01% accuracy criterion that we advocate here is by no means definitive. We 
have seen, however, that even an order of magnitude relaxation in this criterion leads to 
obvious effects on the large-scale field. Hence, it is the least that one should require for a 
transition simulation. 


A Re = 5000 Channel Transition 


Gilbert and Kleiser (1986) have simulated channel flow transition at Re = 5000. Their 
calculation was performed on a 64 x 128 x 128 grid and they were able to maintain an 8 order 
of magnitude reduction in the spectra until t = 134. The results of the previous section 
indicate that greater resolution was required at the lower Reynolds number of 1500. This 
seems counter-intuitive since one expects smaller scales to develop at the higher Reynolds 
number. One reason behind the need for greater resolution for the Re = 1500 case is that 
the initial amplitudes are 3 times larger at the lower Reynolds number. Nonlinear effects 
lead to a more rapid saturation of the higher Fourier-Chebyshev harmonics. Moreover, the 
Chebyshev distribution of grid-points in the normal direction is better tuned to the flow 
field at Re = 5000, which is more concentrated towards the walls. 

These conclusions are supported by our own calculation under the Gilbert and Kleiser 
conditions. The fundamental wavenumbers were a = 1.12 and /? = 2.10. The initial 
amplitudes were = 0.03 and e 3 £> = 0.001, and the initial phase shift 0 = 0°. They used 
the constant flow rate formulation rather than the constant pressure gradient formulation 
employed in the our Re = 1500 calculations. (As discussed by them, the flow settles down 
to a statistically steady turbulent state sooner under the constant flow rate condition.) 
In our repeat of this calculation the constant flow rate formulation was used, and it was 
run on a 64-bit machine. Compared with the case discussed in the previous section, this 
transitional flow takes a lot longer to reach a strongly nonlinear stage since it began from 
much smaller initial amplitudes. 

Figure 8 shows the vertical shear and the streamwise vorticity at t = 128, t=134, 
and t = 138. The resolution spectra at these times are given in Figure 9. There are no 
visible oscillations at t = 128 and t = 134 and the resolution spectra exhibit the 8 order 
of magnitude decay that was found necessary at Re = 1500. This criterion is not satisfied 
at t = 138, particularly in the x-direction, and the oscillations in the top right frame of 
Figure 8 are the consequence. We find that, strange as it may seem at first sight, it is 
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easier to simulate transition at Re = 5000 than at Re = 1500. 

Notice that the curl in the tails of the resolution spectra are much milder here than 
they were for the Re = 1500 case. As documented by Zang (1989), this is one benefit of 
the skew-symmetric form of the Navier-Stokes equations. 


A Re = 1100 Boundary-Layer Transition 


Recently, Zang, Hussaini, and Erlebacher have performed a high resolution simulation 
of boundary-layer transition under the conditions of the experiment by Kovasznay, Ko- 
moda, and Vaseduva (1962). (See Zang and Hussaini (1987) for some preliminary results.) 
The Reynolds number, based on displacement thickness and the free stream velocity, was 
1100, and the fundamental wavenumbers were a = 0.250 and /? = 0.209. The initial am- 
plitudes were €2 d = 0.018 and 630 = 0.008, and the initial phase shift 0 = -36°. The 
simulation was performed under the parallel flow assumption in which the non-parallel ef- 
fects are ignored. This permits much better resolution of the flow at the price of neglecting 
some of the physics of a true, growing boundary layer (see Zang and Hussaini (1987)). 

The calculation was performed on a 144 x 144 x 324 grid and Zang, Hussaini, and 
Erlebacher reproduced numerically the hydrogen-bubble flow visualization results of Hama 
and Nutant (1963). They maintained an 8 order of magnitude reduction in the spectra 
throughout the calculation. Figure 10, taken from that work, displays the vertical shear 
and the mean velocity profiles in the peak plane. The horizontal direction encompasses 
one period of the TS wave and the vertical direction covers 6 displacement thicknesses of 
the undisturbed boundary layer. The vertical shear clearly indicates the primary roll-up of 
the shear layer. (The times in the figure are given in units of the period of the 2-D wave.) 
By Frame 6 a secondary roll-up is just commencing at a height of about 2.5 displacement 
thicknesses. 

Some questions have arisen as to how much resolution is required to reproduce the 
essential features of this flow, at least in the peak plane. This is a particularly important 
question for true spatial simulations because the limitations of current computers mandate 
that resolution in the spanwise direction be sacrificed to furnish the much greater resolution 
needed in the streamwise direction for the spatially-growing boundary layer. To date 
spatial simulations have employed drastic spanwise truncations, retaining only between 4 
and 32 points in z (see Fasel, Rist, and Konzelmann (1987)). With so few points in the 
spanwise direction, truly small-scale spanwise oscillations will not occur - they will be on 
the scale of the spanwise grid, which in these cases is quite coarse. One might hope that 
such simulations will still capture the essential physics at the later stages of transition. 

To investigate this issue, we repeated the Zang, Hussaini and Erlebacher parallel 
boundary-layer simulation from t = 3 onwards using a collocation grid of 96 x 96 x 12. 
To increase the accuracy of this drastically truncated simulation, de-aliasing was applied 
in the x and z directions using the 2/3-rule (see Canuto et al (1988, Chapter 3)). The 
results are shown in Figure 11. The contrast with the well-resolved flow is substantial. 
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The detached shear layer is much less intense, it stays much lower in the boundary layer, 
it lags behind its true streamwise position, and its roll-up is quite mild. 

These plots permit some inferences about what part of the physics of transition is 
mis-represented by the coarse-grid simulation. A principal effect in this type of transition 
(see Zang and Krist (1989)) is the pinching of the lambda vortex towards the peak plane 
and its lifting away from the wall and towards the free stream. Much of this behavior is 
due to the self-induction of the lambda vortex. The intensity of the self-induction is tied 
to the strength of the vorticity. On a coarse grid, only a relatively weak lambda vortex 
is supportable. Hence, the evolution of the lambda vortex is retarded, both in the rate 
at which it pulls away from the wall and the rate at which it moves in the downstream 
direction. Since the detached shear layer is tied to the top of the lambda vortex, its 
intensity, location, and speed is affected in like manner. Figures 10 and 11 are but a 
dramatic illustration of this effect. 

In this simulation even though the resolution in z was woefully inadequate, the resolu- 
tion in x and y was kept sufficient to maintain 0.01% accuracy in those directions. The fact 
that fewer grid-points were required here reflects the phenomenon that since only large 
scales are supported in z and the scales in the three coordinate directions are inter-related, 
fewer small-scale features appear in the other directions. The unfortunate fact is that the 
more grid-points that are used in the spanwise direction, the more that are also needed in 
the other two directions. 

These two simulations beg the question of just how much spanwise resolution is enough. 
Figures 12-14, which contain the results of simulations employing 18, 24, and 48 spanwise 
grid-points, respectively, provide further information. The results clearly improve, i.e., 
tend more towards the well-resolved case, as the number of spanwise grid-points increases. 
All of these intermediate cases exhibit the primary roll-up of the shear layer. However, 
these also display a premature secondary roll-up. 

The principal flaw in them appears to be the occurrence of regions of intense reverse 
spanwise vorticity. The mean wall shear is 0.57, the peak value in the detached shear 
layer of the well-resolved case at Frame 6 is 2.4, and the minimum value is -0.2. The color 
bar was cut off at du/dy < —0.5. Values less than -0.5 appear as white regions. These 
occur below and just upstream of the secondary roll-up in the shear layer. They are as 
strong as -1.0, -1.5, and -1.5 in the respective Frame 6 portions of Figures 12, 13, and 14. 
(Although these figures display the vertical shear, they correspond closely to the (negative 
of the) spanwise vorticity.) Thus, all the intermediate cases are corrupted by the presense 
of incorrect regions of reversed spanwise vorticity. In the case N x = 48 these is additional 
corruption in the near wall region. This reversed vorticity has as yet undetermined effects 
on the physics of the flow but appears related to an exaggeration of the shear-layer roll-up. 


Discussion 


The emphasis in this paper has been on inferring what degree of resolution is needed 
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to capture reliably the instantaneous structure of a transitional flow. This requirement 
stems from the desire to use numerical simulations to clarify the details of the final stages of 
transition and to assist in the development of reliable models for this important flow regime. 
These are more challenging goals numerically than simulating the low-order statistics of 
a turbulent steady state. Direct numerical simulations of turbulent channel flow have 
typically settled for a 3 or 4 order of magnitude reduction in the one-dimensional energy 
spectra (see, for example, Kim, Moin, and Moser(l987)). Turbulence simulations have 
even produced reasonable mean flow and velocity fluctuation profiles in channel flow on 
an 8 x 32 x 8 grid (Gilbert and Kleiser (1986)). These results, however, were obtained 
with the assistance of a temporal average of the computed results over a substantial time 
interval. 

This option is not available for transition simulations. Some alternatives to the straight- 
forward simulations presented here are the incorporation of suitable models of subgrid-scale 
transition structures and the development of efficient self-adaptive numerical algorithms 
that automatically concentrate the grid-points in regions of small-scale structure. Neither 
of these developments, however, is an easy task. 


Conclusion 


Viable short-cuts for reliable direct simulations of transition to turbulence are hard to 
come by. 
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Captions 

Figure 1. Schematic of the channel flow problem. The streamwise coordinate is denoted 
by x, the coordinate normal to the walls by y, and the spanwise coordinate by z. 

Figure 2. Vertical shear (left) and streamwise vorticity (right) contours in the lower half 
of the Re = 1500 channel at t = 15. The contour interval is 0.40. Dashed lines indicate 
negative contours. 

Figure 3. Vertical shear (left) and streamwise vorticity (right) contours in the lower half 
of the Re = 1500 channel at t = 18.75. 

Figure 4. Resolution spectra in each of the coordinate directions for Re = 1500 channel 
flow simulations. 

Figure 5. Vertical shear (left) and streamwise vorticity (right) contours in the lower half 
of the Re = 1500 channel at t = 18.75 for a 64 s FFDM simulation. 

Figure 6. Vertical shear (left) and streamwise vorticity (right) contours in the lower half 
of the Re = 1500 channel at t = 22.5. 

Figure 7. Vertical shear (left) and streamwise vorticity (right) contours in the lower half 
of the Re = 1500 channel at t = 27.5. 

Figure 8. Vertical shear (top) and streamwise vorticity (bottom) contours in the lower half 
of the Re = 5000 channel. The contour interval is 0.40. 

Figure 9. Resolution spectra in each of the coordinate directions for the Re = 5000 channel 
flow simulation. 

Figure 10. Vertical shear contours and mean velocity profiles for the Re = 1100 boundary- 
layer simulation on a 144 x 144 x 324 grid. 

Figure 11. Vertical shear contours and mean velocity profiles for the Re = 1100 boundary- 
layer simulation on a 96 x 96 x 12 grid. 

Figure 12. Vertical shear contours and mean velocity profiles for the Re = 1100 boundary- 
layer simulation on a 192 x 128 x 18 grid. 

Figure 13. Vertical shear contours and mean velocity profiles for the Re = 1100 boundary- 
layer simulation on a 192 x 96 x 24 grid. 

Figure 14. Vertical shear contours and mean velocity profiles for the Re = 1100 boundary- 
layer simulation on a 256 x 128 x 48 grid. 
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